package scu.maqiang.elasticity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q42;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.numeric.*;

public class Elasticity2DQ42 {

	public static void main(String[] args) {
//	    use SMatrix_m
//	    use Mesh2DQ4_m
//	    use FES2DQ42_m
//	    use IterSSolver_m
//	    implicit none
//	    
//	    type(mesh2DQ4) :: mesh2D
//	    type(FES2Q42) :: fs
//	    type(SMatrix) :: A
//	    type(IterSSolver) :: solver
//	    real*8, pointer :: x(:), RHS(:), u(:), v(:)
//	    real*8, pointer :: ex(:), ey(:), exy(:)
//	    real*8, pointer :: tx(:), ty(:), txy(:)
//	    integer:: label(1)
//	    integer :: verbosity = 1
//	    character*80 :: outputFile

	    Mesh2Q4 mesh2D = new Mesh2Q4().square2D(80, 10);
	    mesh2D.scale(8.0, 1.0);
	    FES2Q42 fs = new FES2Q42(mesh2D);
	    SRMatrix A = new SRMatrix(fs.GetNdof());
	    double[] ENu = {1.0e6, 0.3};
	    fs.assembleStiff(ENu, BVPType.PLANESTRAIN, A);
	    
	    double[] RHS = new double[fs.GetNdof()];
	    double[] force = {100.0, 0.0};
	    fs.assembleFlux(force, BVPType.COMMON, RHS, 2);

	    System.out.println("Before Apply BC");
		System.out.println("The Size Of SparseMatrix " + A.getNNZ());
	    
		fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 4);
	    
		System.out.println("After Apply BC");
		System.out.println("The Size Of SparseMatrix " + A.getNNZ());
		double[] x = new double[fs.GetNdof()];
		
//		NewIterSSolver solver = new NewIterSSolver(A);
//		solver.PCGSSOR(RHS, x, 1.5, 1);
		DirectSSolver solver = new DirectSSolver(A);
		x = solver.CholeskySolve(RHS, 1);
		double[][] disp = new double[2][fs.GetNdof() / 2];
		fs.extractComponent(x, disp);
		
	    System.out.println("Displacement:");
		System.out.println("u min : " + MVO.min(disp[0]) + " max: " + MVO.max(disp[0]));
		System.out.println("v min : " + MVO.min(disp[1]) + " max: " + MVO.max(disp[1]));
		String outputFile = "solutionPlainStress2DR.dat";
		mesh2D.toTecplot(outputFile, disp);
		String deformFile = "Q4Deform.dat";
		mesh2D.toTecplotDeform(deformFile, disp, 1.0e4);
		
//	    
//	    call fs.computeStrain(u, v, ex, ey, exy)
//
//	    write(*, *) "Strain: "
//	    write(*, *) "ex min : ", minval(ex), "max: ", maxval(ex)
//	    write(*, *) "ey min : ", minval(ey), "max: ", maxval(ey)
//	    write(*, *) "exy min : ", minval(exy), "max: ", maxval(exy)
//
//	    !call fs.computeStressQ42(ex, ey, exy, funcE, funcNu, tx, ty, txy, 0)
//	    !write(*, *) "Stress: "
//	    !write(*, *) "tx min : ", minval(tx), "max: ", maxval(tx)
//	    !write(*, *) "ty min : ", minval(ty), "max: ", maxval(ty)
//	    !write(*, *) "txy min : ", minval(txy), "max: ", maxval(txy)
//
//	    outputFile = "solutionPlainStress2DR.dat"
//	!    call exportMeshSolution(outputFile, mesh2D, u, v)
//	!    outputFile = "solutionPlainStress2DRStrain.dat"
//	!    call exportMeshSolution(outputFile, mesh2D, ex, ey, exy)
//	!    outputFile = "solutionPlainStress2DRStress.dat"
//	!    call exportMeshSolution(outputFile, mesh2D, tx, ty, txy)
	}

}
